A Choropleth map is a thematic map in which areas are colored considering the quantitative measurement of a variable, such as population density or gross domestic product.
There are a lot of libraries in R (and of course in Python too) that perform it efficiently. However, I feel that ggplot2 offers the fastest and most friendly way. Don’t you believe me? Look at this.
geom <- st_read(system.file("shape/nc.shp", package="sf"))
ggplot(geom) +
geom_sf(aes(fill = AREA)) +
ggtitle('My first choropleth map') #so intuitively, right?
ggplot2 has a lot of tricks and well-written documentation, however, it does not support dynamics and interactive visualization natively. To accomplish this task, one option is changing ggplot by plotly, but don’t worry ggplot2 and plotly working very well together and are similar, because both are based on the grammar of graphics.
geom <- st_read(system.file("shape/nc.shp", package="sf"),quiet = TRUE)
p <- ggplot(geom) +
geom_sf(aes(fill = AREA)) +
ggtitle('My first choropleth map')
ggplotly(p)
In this post, I aim to show you how to construct a beautiful 3D interactive choropleth map in just a few minutes using R.
Plotly or Plot.ly, is a technical computing company that provides scientific graphing libraries for different languages. Theses libraries use as a basis plotly.js a Javascript library building on top of d3.js and stack.gl. In the case of R, the plotly library take advantage of the htmlwidget framework (see plotly::as_widget) to create a “connection” with plotly.js. You can install plotly in R as follow:
devtools::install_github("ropensci/plotly")
In a nutshell, there are four ways to create geospatial visualization from a sf object.
plot_ly( ): Is the main function that pass R objects (list) to plotly.js.
ggplotly( ): This function converts a ggplot2::ggplot() object to a plotly object. Does not support orthographic (3D) projections.
plot_mapbox( ): It is similar to use plot_ly(type = 'scattermapbox'), allow customization of basemaps.
plot_geo( ): It is similar to use plot_ly(type = 'scattergeo').
If you are trying to visualize a plot with a lot of graphical elements to render, I highly recommend use plot_ly instead of ggplotly. plot_ly is the most stable procedure to leverage your plot to the plotly.js API. For a large explanation about the strengths and weaknesses of referred functions read this section of the plotly ebook.
library("rnaturalearth")
library("classInt")
library("stringr")
library("viridisLite")
library("tidyverse")
library("sf")
library("smoothr")
library("plotly")
For download the world population dataset, I used the rnaturalearth::ne_download function. It return a sf object in memory if load = TRUE and returnclass = "sf".
worldmap <- ne_download(scale = 110,
type = "countries",
category = "cultural",
destdir = tempdir(),
load = TRUE,
returnclass = "sf")
The sf object worldmap has 177 rows (countries) and 95 columns (variables), we selects the POP_EST (population) and NAME (countries names) columns.
worldmap_pop <- worldmap %>%
select('POP_EST', 'NAME') %>% # country population
'colnames<-'(c('pop','name','geometry')) %>%
mutate(pop = round(as.numeric(pop)/1000000,2)) %>%
arrange(pop)
# Filling polygons with holes. (Sud-Africa case)
sudafrica <- worldmap_pop[153,] %>%
fill_holes(10^12) %>%
st_cast('MULTIPOLYGON')
worldmap_pop[153,] = sudafrica
intervals <- classIntervals(worldmap_pop$pop,style = 'kmeans',n = 10)$brks
intervals_f <-cut(worldmap_pop$pop,intervals)
intervals_f <- factor(intervals_f,rev(levels(intervals_f)))
lvls <- levels(intervals_f) %>%
str_replace_all("\\(|\\]","") %>%
str_replace_all(","," - ")
lvls[c(1,length(lvls))] <- c('804 >', '< 4.08')
levels(intervals_f) <- lvls
worldmap_pop$interval <- intervals_f
text <- sprintf('Country: %s \nPopulation: %.2f',worldmap_pop$name,worldmap_pop$pop)
worldmap_pop$text =text
# geospatial parameters --------------------------------------------------
geo <- list(showland = FALSE,
showlakes = FALSE,
showcountries = TRUE,
showocean = TRUE,
countrywidth = 0.5,
landcolor = toRGB("grey90"),
lakecolor = toRGB("white"),
oceancolor = toRGB("#e5f3fc"),
projection = list(type = 'orthographic',
rotation = list(lon = -100,
lat = 40,
roll = 0)),
lonaxis = list(showgrid = TRUE,
gridcolor = toRGB("gray80"),
gridwidth = 0.5),
lataxis = list(showgrid = TRUE,
gridcolor = toRGB("gray80"),
gridwidth = 0.5))
# Legend parameters -------------------------------------------------------
legend_param <- list(x = 1.1,
y = 0.82,
font = list(family = "sans-serif",
size = 18,
color = "black"),
bgcolor = "#E2E2E2",
bordercolor = "#FFFFFF",
borderwidth = 2)
# Image title parameters -------------------------------------------------
legend_title <- list(yref = "paper",
xref = "paper",
y = 0.98,
x = 1.4,
text = "<b>World population \n in millions</b>",
font = list(color = "black",
family = "sans serif",
size = 22),
showarrow = FALSE)
p <- plot_ly(data = worldmap_pop,
text = ~text,
color = ~interval,
colors = viridis(10,direction = -1), # the number of colors need to be according to intervals (classInt::classIntervals)
type = "scattergeo",
alpha = 1,
stroke = I("black"), # boundary color line
hoverinfo = "text", # Display just the ~text
span = I(0.5)) %>% # thickness of boundary line
hide_colorbar() %>%
layout(showlegend = TRUE,
geo = geo,
legend = legend_param,
annotations = legend_title) %>%
config(displayModeBar = FALSE) #do not display plotly toolkit.
p
Plotly is an incredible package for making interactive visualization, however, the enormous amount of graphical elements could be overwhelming. The solution is to create Code snippets.
Code snippets are text macros that are used for quickly inserting common pieces of code. In Rstudio, you can edit the built-in snippet definitions and even add snippets of your own via the Edit Snippets button in:
Tools -> Global Options -> Code -> Edit Snippet
My snippet for making 3D choropleths maps (worldmap) are defined below:
snippet worldmap
# geospatial parameters --------------------------------------------------
geo <- list(showland = FALSE,
showlakes = FALSE,
showcountries = TRUE,
showocean = TRUE,
countrywidth = 0.5,
landcolor = toRGB("grey90"),
lakecolor = toRGB("white"),
oceancolor = toRGB("#e5f3fc"),
projection = list(type = 'orthographic',
rotation = list(lon = -100,
lat = 40,
roll = 0)),
lonaxis = list(showgrid = TRUE,
gridcolor = toRGB("gray80"),
gridwidth = 0.5),
lataxis = list(showgrid = TRUE,
gridcolor = toRGB("gray80"),
gridwidth = 0.5))
# Legend parameters -------------------------------------------------------
legend_param <- list(x = 1.1,
y = 0.82,
font = list(family = "sans-serif",
size = 18,
color = "black"),
bgcolor = "#E2E2E2",
bordercolor = "#FFFFFF",
borderwidth = 2)
# Image title parameters -------------------------------------------------
legend_title <- list(yref = "paper",
xref = "paper",
y = 0.98,
x = 1.4,
text = "<b>${1:Title}</b>",
font = list(color = "black",
family = "sans serif",
size = 22),
showarrow = FALSE)
p <- plot_ly(data = ${2:yourdata},
text = ~${3:text},
color = ~${4:color},
colors = viridis(${5:N},direction = -1),
type = "scattergeo",
alpha = 1,
stroke = I("black"),
hoverinfo = "text",
span = I(0.5)) %>%
hide_colorbar() %>%
layout(showlegend = TRUE,
geo = geo,
legend = legend_param,
annotations = legend_title) %>%
config(displayModeBar = FALSE)
p